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Abstract 

This  paper  presents  a  series  of  tests  that  have  been  conducted  to  validate  the  collision 
models  of  the  H aw  kDSMC  implementation.  These  tests  cover  the  hard  sphere  (HS),  vari¬ 
able  hard  sphere  (VHS),  and  variable  soft  sphere  (VSS)  collision  models,  and  the  Larsen- 
Borgnakkc  energy  exchange  model,  for  non-reacting  gas  mixtures.  The  results  show  that 
//mvTobtains  the  correct  results  in  all  cases  under  consideration. 


1  Introduction 

Recent  advances  in  microprocessor  performance  have  been  driven  primarily  by  improvements 
in  manufacturing  technology.  New  processes  and  equipment  have  paved  the  way  for  smaller 
feature  sizes  and  larger  die  sizes.  These  have  in  turn  enabled  the  production  of  microprocessors 
with  more  transistors,  operating  at  lower  voltages  and  higher  clock  rates.  One  of  the  key  pieces 
of  equipment  in  microelectronics  manufacturing  is  the  plasma  reactor ,  used  in  30  to  40  percent 
of  the  processing  steps.  Plasma  reactors  use  energetic  rarefied  gases,  plasmas,  to  remove  parti¬ 
cles  from,  and  deposit  particles  on,  silicon  wafers.  Improving  the  design  of  these  reactors,  and 
the  processes  that  they  are  used  for,  will  enable  the  microelectronics  industry  to  make  smaller, 
cheaper,  and  faster  microprocessors. 

Design  and  optimization  of  plasma  reactors  has  so  far  been  largely  empirical.  Experiments 
have  been  conducted  to  improve  process  configurations,  but  because  of  the  high  equipment  and 
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Figure  1 :  A  Plasma  Reactor 


operating  costs  of  plasma  reactors,  detailed  parametric  studies  have  been  economically  imprac¬ 
tical.  Computer-based  simulation  of  the  plasma  flow  inside  a  reactor  will  allow  manufacturers 
to  evaluate  the  viability  of  different  reactor  designs  before  they  are  implemented.  Once  a  reac¬ 
tor  has  been  installed,  simulation  results  will  also  be  useful  for  studying  the  effects  of  different 
operating  conditions,  thereby  optimizing  processing  stages. 

Figure  1  shows  a  typical  plasma  reactor,  the  GEC  Reference  Cell,  and  Figure  2  schemati¬ 
cally  depicts  its  operation.  A  silicon  wafer  is  attached  to  an  electrode,  and  plasma  fills  the  space 
between  the  wafer  and  another  electrode.  Gas  flows  in,  reactions  take  place  within  the  gas  and 
on  the  surface  of  the  wafer,  and  the  products  of  these  reactions  are  pumped  out  of  the  reactor. 
Electromagnetic  fields,  applied  through  the  electrodes,  add  energy  to  the  system. 

Advanced  plasma  simulation  capabilities  will  be  directly  applicable  to  problems  in  the  mi¬ 
croelectronics  industry  and  can  therefore  have  direct  bearing  on  industrial  competitiveness.  Sim¬ 
ulations  will  be  useful  for  studying  process  optimization,  compact  model  development,  equip¬ 
ment  evaluation,  process  control,  and  technology  feasibility.  Efficient  modeling  will  reduce  the 
time  and  cost  of  microelectronics  development,  and  therefore  help  improve  the  quality  of  the 
next  generation  of  microprocessors. 
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Figure  2:  Reactor  Schematic 


2  Direct  Simulation  Monte  Carlo 


Plasma  flow  is  in  the  transition  regime :  the  mean  free  path  of  particles  is  too  large  for  traditional 
continuum  CFD  methods  to  be  applicable,  but  because  collisions  are  important,  free-molecular 
simulations  are  not  appropriate,  either.  The  Direct  Simulation  Monte  Carlo  (DSMC)  method 
(See,  for  example,  [Bird94].)  is  a  widely  used  numerical  approach  for  solving  rarefied  gas  dy¬ 
namics  problems  in  the  transitional  regime.  It  simulates  individual  particles  as  they  move  through 
space  and  collide  with  solid  objects  and  other  particles.  Macroscopic  properties,  such  as  density 
and  temperature,  can  be  computed  by  appropriate  averaging  of  particle  masses,  positions,  and 
velocities.  Surface  properties  are  calculated  from  the  momentum  and  energy  exchanges  during 
collisions  with  surfaces. 

In  a  DSMC  simulation,  a  physical  region  of  interest  is  decomposed  into  a  number  of  cells. 
The  cells  are  initially  filled  with  simulation  particles  according  to  density,  temperature,  and  ve¬ 
locity  specifications.  The  simulation  takes  discrete  steps  in  time,  during  which  these  particles  are 
allowed  to  move  throughout  the  domain  and  collide  with  other  particles.  Figure  3  shows  several 
cells  (in  two  dimensions)  and  some  of  the  possible  operations  that  can  take  place  on  a  particle 
during  a  timestep.  For  example,  particles  may  flow  into  the  domain  through  injection  cells,  or 
out  of  the  domain  through  exhaust  cells.  They  may  also  collide  with,  or  become  embedded  in, 
solid  objects  in  the  domain. 
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Particle  Leaving  Domain 


Particle  Moving  Between  Cells 


Particle  Emitted  by  Surface 


Particle  Absorbed  by  Surface 


Figure  3:  DSMC  Cells  and  Particles 


The  DSMC  method  simulates  the  Boltzmann  equation  by  decoupling  particle  movement  and 
collisions.  The  two  main  components  of  any  DSMC  implementation  are  therefore  the  transport 
and  collision  algorithms.  Particle  transport  is  conceptually  straightforward,  and  is  discussed  in 
detail  in  [Rieffel95].  Both  the  power  and  computational  cost  of  the  DSMC  method  are  due  to 
the  collision  process,  and  are  therefore  the  subject  of  this  paper. 

The  main  goal  of  this  work  is  to  validate  the  collision  models  that  have  been  implemented 
in  Hawk ,  a  DSMC  application  designed  for  the  simulation  of  plasma  reactors.  General  charac¬ 
teristics  of  Hawk,  including  methods  for  computing  macroscopic  parameters,  grid  issues,  and 
the  parallel  implementation,  are  discussed  in  [Rieffel95].  The  main  purpose  of  this  paper  is  to 
present  the  collision  models  that  are  implemented  in  Hawk,  and  to  demonstrate  their  validation. 

For  the  purposes  of  validation,  it  has  been  useful  to  compare  the  results  obtained  with  Hawkto 
those  obtained  with  another  DSMC  implementation,  SMILE.  SMILE  is  a  Computational  tool 
for  solving  problems  of  rarefied  gas  aerodynamics  [Ivanov92],  created  at  the  Institute  of  The¬ 
oretical  and  Applied  Mechanics  (Novosibirsk,  Russia).  It  is  based  on  the  majorant  principle  of 
construction  and  substantiation  of  numerical  schemes  for  the  DSMC  method.  The  coupling  of 
’’cell”  and  ’’free  cell”  schemes  [Ivanov94]  provides  the  required  spatial  resolution  throughout 
the  whole  flow  field,  including  regions  with  strong  gradients.  The  preprocessing  subsystem  of 
this  tool  is  used  for  defining  the  geometric  model  of  a  space  vehicle  and  for  specifying  bound¬ 
ary  and  initial  conditions.  The  results  of  a  computation,  both  flow  fields  and  distributed  surface 
characteristics,  can  be  analyzed  with  the  postprocessing  system. 

The  following  sections  describe  the  collision  models  that  have  been  implemented  in  Hawk, 
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and  present  results  of  a  series  of  validation  tests,  comparing  against  analytical  results  as  well  as 
the  results  of  SMILE. 
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3  DSMC  Collision  Overview 


Generally,  the  global  problem  of  creating  collision  models  for  DSMC  method  may  be  divided 
into  several  subproblems,  related  to  creation  of 

-  intermolecular  potential  models 

-  models  of  rotational  degrees  of  freedom  of  molecules 

-  models  of  vibrational  degrees  of  freedom  of  molecules 

-  models  of  chemical  reactions. 

The  basic  models  for  all  four  groups  used  in  DSMC  method  are  presented  below. 


3.1  The  Lennard-Jones  Potential 


A  number  of  models  to  describe  intermolecular  potential  is  suggested  in  references.  It  is  usually 
assumed  that  the  potential  function  includes  molecular  attraction  at  large  distances  and  repulsion 
at  small  ones.  The  most  famous  and  mathematically  convenient  potential,  taking  into  account 
molecular  attraction  and  repulsion,  is  the  Lennard-Jones  potential, 

<f>(r)  —  4e 


a 


12 


a 


where  a  is  the  distance  at  which  the  potential  function  changes  its  sign,  e  is  the  minimum  poten¬ 
tial  value,  and  r  is  the  distance  between  the  particle  centers.  The  deflection  angle,  x,  is, 

fWi 

I  o 


X  =  7T  -  2 

where  W\  >  0  is  the  root  of  equation, 


1  -  W2  - 


\mR92 , 


dW, 


1-W2-  =  0, 

2  mR9 

W  =  b/r.  and  b  is  the  distance  of  closest  approach  of  the  undisturbed  molecular  trajectories  in 
the  center  of  mass  frame  of  reference. 

In  dimensionless  variables,  x  may  be  rewritten  as 
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_  mR9 
2e 


where  W  =  £,  W0  =  6*  (£)  12  ,  IV,'  =  6*  (g)  \  6*  =  r*  =  f,  and  g *  = 

In  statistical  simulation,  the  dependence  x  on  two  parameters  Wo  and  IV,'  necessitates  the  use 
of  interpolation  tables  for  computing  the  maximum  impact  parameters  and  deflection  angles.  To 
find  the  deflection  angle  for  a  pair  of  colliding  particles  in  simulation,  it  is  necessary  to  determine 
the  value  b *max  from  the  interpolation  tables  at  the  specified  deflection  angle  cut-off,  sample  6* 
and,  using  the  tables  again,  find  y.  Because  of  the  high  computational  cost  of  this  process,  as 
well  as  difficulties  in  determining  the  mean  free  path,  the  Lennard-Jones  potential  is  typically 
only  applied  to  homogeneous  relaxation  cases  and  simple  one-dimensional  problems. 


3.2  The  Inverse  Power  Law  Potential 

A  more  suitable  model  for  high-temperature  flows,  where  intermolecular  repulsion  prevails,  is 
the  Inverse  Power  Law  (IPL)  model, 


<f>(r 


(//  —  ’ 

where  g  is  the  power  determining  the  ’’hardness”  of  particles.  The  deflection  angle  is, 


X  =  *  ~  2 


fWi 


1  -  W2  - 


g  -  1  \WoJ 


wy-1 


dW, 


where  W0  —  b  1,-1 ,  and  \\\  >  0  is  the  root  of  equation  1  —  W2  —  )  V  =0.  Note 

that  x  is  a  function  of  one  parameter,  Wo,  and  may  be  found  from  various  approximate  analytical 
expressions  [Nanbu81]. 

The  total,  viscosity,  and  diffusion  cross-sections  for  IPL  model  are, 


<7j  =  j  crdQ  =  oo 


cr, 


1  -  cos2x)vdn  =  2t tA2(v)  j 


<jd  —  /  ( 1  —  cosx)o'd9t  —  2'kAi(v) 


v- 1 


m.RK 

y  92  J 

To  use  the  IPL  model  in  a  DSMC  simulation,  it  is  necessary  to  limit  crj .  This  can  be  achieved 
by  specifying  either  a  deflection  angle  cut-off  or  a  maximum  impact  parameter  value.  In  most 
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cases,  the  latter  is  chosen.  Since  \  is  a  function  of  W0 ,  the  deflection  angle  cut-off  may  be  applied 
through  the  specification  of  a  maximum  value  Wo,m  of  Wq.  This  yields, 


=  7T  lly 


rj-1 


0  ,m  \  9 

\m.Rgz 


At  Tj  —  oo,  the  IPL  model  is  equivalent  to  the  hard  sphere  model  for, 


<tj  =  7 rd2,  b  —  dcos{x/2 ), 


where  d  is  the  molecular  diameter. 
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4  The  Hard  Sphere  Collision  Model 

Earlier  versions  of  this  implementation  have  been  validated  by  comparison  with  published  re¬ 
sults  [Zhong95]  for  heat-transfer  problems  [Rieffel95].  This  section  presents  several  further 
tests,  using  the  Hard  Sphere  (HS)  collision  model  for  gas  mixtures  and  more  complicated  ge¬ 
ometries. 

4.1  Box  Tests 

The  first  series  of  tests  considered  here  was  designed  to  verify  that  all  of  the  basic  DSMC  oper¬ 
ations  were  correctly  implemented,  including  transport,  inflow,  accomodation,  and  hard  sphere 
(HS)  collisions.  This  series  was  performed  in  a  12-cell  uniform  cubic  grid,  with  all  walls  spec¬ 
ular. 

The  first  test  considered  Ar  and  Ar*,  two  identical  species.  It  was  verified  that  the  collision 
frequencies  were  correct.  Macroparameters  were  observed  to  be  identical  for  both  species. 

The  second  test  considered  Ar  and  He  in  equal  concentrations,  with  density  n  =  1.83  x  1019 
particles/m3,  and  both  species  at  a  temperature  of  300K.  The  system  remained  in  equilibrium, 
and  all  collision  counts  were  equal  to  theoretical  values. 

The  next  test  considered  a  mixture  of  Argon  and  Helium  at  different  temperatures.  The  sys¬ 
tem  reached  equilibrium  of  Tat=  Tjje  =  300K,  the  collision  numbers  were  correct,  and  the  system 
reached  equilibrium  at  the  same  rate  as  SMILE. 

The  next  test  considered  different  concentrations  of  Ar  and  He,  [Ar]=0.9,  [He]=0.1.  Tat  = 
100K  7  j ie  =  500K.  Verified  that  the  final  temperature  for  both  species  was  0.9*100  +  0.1*500  = 
140K. 

The  last  box  case  considered  the  opposite  concentrations,  [Ar]=0.1,  [He]=0.9,  with  TAt  = 
100K  and  J'u,  =  500K.  The  final  temperature  for  both  species  was  the  correct  value,  0.1  x  100  H- 
0.9  x  500  =  460 Ji. 

4.2  Cylinder  Test 

The  final  HS  test  considers  Mach-4  flow  past  a  cylinder.  The  cylinder  has  a  radius  of  0.05m, 
and  the  grid  extends  0. 15m  in  front  of  the  center  of  the  cylinder,  0.35m  behind  the  center  of  the 
cylinder,  and  0.2m  to  the  side  of  the  cylinder.  The  freestream  density  is  1.83  x  1021  particles  per 
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cubic  meter,  the  freestream  speed  is  526.86  m/s,  and  the  freestream  temperature  is  50  K.  The 
cylinder  itself  is  fully  accomodating  at  300K.  The  inflow  is  composed  of  50%  Ar  and  50%  He. 

Because  the  symmetric  nature  of  the  problem,  only  the  upper  half  of  the  cylinder  was  simu¬ 
lated.  A  thin  (0.01m)  3-D  grid  was  used  for  Hawk ,  while  a  2D  grid  was  used  for  SMILE.  Argon 
density  plots  for  both  codes  are  shown  in  Figures  4  and  5.  The  results  of  the  simulations  agree 
very  well. 
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Figure  4:  HawkAr  Density 


Figure  5:  SMILE  Ar  Density 
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5  The  Variable  Hard  Sphere  Collision  Model 


The  Variable  Hard  Sphere  (VHS)  model,  developed  by  Bird  in  1981,  uses  cross  sections  that  are 
functions  of  relative  velocity,  but  with  hard  sphere  (HS)  scattering  angles.  Collision,  viscosity, 
and  diffusion  cross  sections  therefore  take  the  form, 


=  7T d2  =  c(mRg2  / 2)  " 

2  72 

a.,  —  —7 rd 

»  3 

(Tjj  —  7T  (l2. 


where  c  and  a  are  parameters  of  the  model.  As  in  the  HS  model,  the  scattering  angle  is  given 

by, 

b  =  dcos(x/2). 


Equating  for  the  IPL  and  VHS  models  yields, 

2 

a  =  - c  =  37rA2(^)(«:/2)a. 

For  a  gas  at  equilibrium,  uj  may  be  written  as, 


&T  —  &ref 


TTBf\ 

~f) 


a 
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where  crref  is  the  reference  value  of  the  cross-section  at  the  reference  temperature  Trej. 

The  VHS  model  is  currently  the  most  popular  for  DSMC  simulations,  because  of  its  sim¬ 
plicity  and  its  good  approximation  to  real  intermolecular  potentials.  Note  also  that  VHS  colli¬ 
sions  are  much  more  ’’efficient”  than  IPL  collisions.  For  example,  relaxation  for  Maxwellian 
molecules  under  the  IPL  model,  with  the  deflection  angle  cut-off  1°,  is  about  six  times  slower 
than  under  the  VHS  model  with  the  same  viscosity-temperature  dependence. 

In  some  cases,  the  parameters  a,  arej  ,  and  Trej  may  be  known  for  collisions  between  par¬ 
ticles  of  a  species  A,  and  similarly  for  a  species  B,  but  not  for  collisions  between  species  A  and 
B.  In  this  case,  they  can  be  approximated  according  to  the  formulae, 


(XAB 


OLA  +  aB 


&AB 


■s/va  +  \/^B^ 


Tab  — 


2 

T A  +  Tb 
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Figure  6:  VHS  50-50 


5.1  Box  Tests 

Initial  validation  of  the  VHS  model  consisted  of  running  a  uniform  specular  box  and  comparing 
the  number  of  collisions  between  Hawkand  SMILE.  The  results  were  in  excellent  agreement. 
The  next  series  of  tests  studied  the  temperature  relaxation  rate  in  the  specular  box,  with  mixtures 
of  Argon  and  Helium  at  different  temperatures  and  in  different  concentrations. 

The  first  of  these  tests  was  a  50%-50%  mixture  of  Argon  and  Helium,  with  Argon  at  100 
K  and  Helium  at  500K.  As  would  be  expected,  the  final  temperature  of  the  mixture  was  300K. 
Figure  6  shows  that  SMILE  and  Hawkboth  reached  the  correct  final  temperature  at  the  same  rate. 

The  second  test  was  with  90%  Argon  at  100K  and  10%  Helium  at  500K.  The  number  of 
collisions,  and  convergence  rates,  as  shown  in  Figure  7,  both  agreed. 

5.2  Cylinder  Test 

The  final  test  of  the  VHS  model  was  the  cylinder  problem,  as  described  in  Section  4.2.  For  this, 
we  used  a  mixture  of  50%  Argon  and  50%  Helium,  with  a=0.5.  Density  plots  for  the  two  codes 
are  shown  in  Figures  8  and  9,  and  these  agree  very  well.  The  density  along  the  streamline  was 
also  compared  and  found  to  be  the  same  for  the  two  codes. 
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Figure  7:  Ar-He  90-10  VHS 
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Figure  8:  Hawk  Ax  Density 


Figure  9:  SMILE  Ar  Density 
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6  The  Variable  Soft  Sphere  Collision  Model 


The  main  disadvantage  of  the  VHS  model  arises  when  equating  viscosity  coefficients  cq,  for  the 
VHS  and  IPL  models.  The  diffusion  cross-sections  for  the  two  models,  and  hence,  the  diffusion 
coefficients,  coincide  only  for  the  hard  sphere  model.  The  variable  soft  sphere  (VSS)  model 
developed  by  Koura  (1991)  has  no  such  drawback.  In  this  model,  the  primary  equations  are, 


<rT  =  7T d2  =  c'(^mRg2)  a 


6ip 


O  +  IXV’  +  2)  3 
2 


-7 Td2 


CTD 


-7 Td2 


i>  + 1 

b  —  dcos^(x/2). 


Equating  the  VSS  and  IPL  parameters,  we  obtain, 
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a  — 


rj  +  r 


i>  = 


'Mv) 

AM 


-i 


c  —  3TrA2(r])(K,/2y 


,(^  +  l)(V>  +  2) 

6  ip 


In  a  collision,  we  define  two  scattering  angles,  x  an(l  e-  X  's  the  angle  between  the  pre¬ 
collision  relative  velocity  and  the  post-collision  relative  velocity,  e  is  the  azimuthal  impact  angle 
measured  between  the  collision  plane  and  some  reference  plane.  For  hard  sphere  collisions  (HS 
and  VHS),  the  scattering  angles  x  and  e  are  both  distributed  uniformly.  In  terms  of  the  molecule 
diameter  d  and  the  impact  parameter  b,  we  can  write 

—  =  cos  — 


d 

For  VSS  colisions,  we  use  a  parameter  §  to  characterize  the  anisotropy  of  the  scattering  angle, 
and  can  therefore  write, 


V 

d 
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COS 


These  parameters  must  be  used  for  computing  the  post-collision  relative  velocity  g'  as  a  func¬ 
tion  of  the  pre-collision  relative  velocity  go  and  the  post-collision  relative  speed  g' .  For  elastic 
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collisions,  g'  is  simply  the  magnitude  of  g,  while  for  inelastic  collisions  it  may  be  different.  The 
procedure  for  computing  g'  is  as  follows. 

First,  the  azimuthal  impact  angle  e  is  computed  using  a  random  number  RL, 

e  =  2ttRi. 


Second,  the  scattering  angle  \  is  computed  using  another  random  number  R2, 


cos  X  =  2 R%  -  1 


sm  x 


yfl 


cosz  x- 


Next,  the  pre-collision  relative  velocity  g0  is  scaled  to  have  the  magnitude  of  the  post-collision 
relative  speed, 


9  =  9 


/  go 


1 00  r 


The  components  of  g'  can  then  be  computed  as 


gx  =  gx  cos  x  +  Jg‘  +  g%  sm  x  sm  e 


g'y  =  gy  cos  x  +  sin  x  ( g' gz  cos  t  -  gxgy  sin  t)  (g2y  +  gfj 


-1/2 


g'z  =  9z  cos  x  -  sin  y  (g'gy  cos  e  +  gxgz  sin  e)  +  gfj 


-1/2 


6.1  Test  1 

The  first  VSS  test  case  was  the  box  problem  with  Helium  and  Argon  parameters  as  shown  in  the 
table  below.  The  number  of  collisions  was  compared  with  SMILE.  The  system  was  shown  to 
stay  at  equilibrium  at  100K. 


Species 

Ar-Ar 

Ar-He 

He-He 

a 

0.31 

0.235 

0.16 

P 

0.714 

0.754 

0.794 

Fraction 

0.5 

0.5 

Temperature 

100 

100 

K 
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T,  K 


Figure  10:  VSS  50-50 


6.2  Test  2 

The  next  test  for  the  VSS  model  was  with  Argon  and  Helium  at  different  initial  temperatures,  as 
shown  in  the  table  below.  The  convergence  rates  of  the  temperatures  were  shown  to  agree  with 
the  results  of  SMILE,  as  shown  in  Figure  10. 


Species 

Ar-Ar 

Ar-He 

He-He 

a 

0.31 

0.235 

0.16 

P 

0.714 

0.754 

0.794 

Fraction 

0.5 

0.5 

Temperature 

100 

500 

K 

6.3  Test  3 

The  final  box  test  for  the  VSS  model  used  Argon  and  Helium  in  different  concentrations,  ac¬ 
cording  to  the  table  below.  Results  were  in  close  agreement  with  those  of  SMILE,  as  shown  in 
Figure  12. 


Parameter 

Ar-Ar 

Ar-He 

He-He 

a 

0.31 

0.235 

0.16 

P 

0.714 

0.754 

0.794 

Fraction 

0.9 

0.1 

Temperature 

100 

500 

K 
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T  (K) 


Figure  11:  VSS  90-10 


Figure  12:  VSS  90-10 
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Figure  13:  Hawk  Ay  Density 


Figure  14:  SMILE  Ar  Density 


6.4  Cylinder 

The  final  VSS  test  considered  the  cylinder  problem,  as  described  in  4.2.  Results  for  Hawkand 
SMILE  are  shown  in  Figures  13  and  14. 
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7  The  Larsen-Borgnakke  Energy  Exchange  Model 


The  Larsen-Borgnakke  model  of  energy  exchange  is  used  to  describe  internal  energy  modes  for 
rotation  and  vibration[Borgnakke75].  Relative  translational  and  internal  post-collision  energies 
of  colliding  particles  are  assumed  to  be  distributed  according  to  equilibrium  distribution  func¬ 
tion. 

The  model  associates  with  each  species  a  number  of  atoms  na  and  a  characteristic  vibrational 
temperature  Qv,  and  with  each  particle,  rotational  energy  Er  and  vibrational  energy  Ev.  Rota¬ 
tional  and  vibrational  energies  are  assumed  to  be  continuous. 

This  implementation  considers  two  types  of  energy  exchange,  translational-rotational  (TR) 
and  translational-rotational- vibrational  (TRV).  Each  collision  has  some  probability  <pj{  of  a  TR 
exchange  and  some  (smaller)probability  (pv  of  a  TRV  exchange. 


7.1  Internal  Degrees  of  Freedom 


It  is  first  necessary  to  characterize  the  number  of  degrees  of  freedom  in  each  of  the  energy  modes. 
Relative  translational  energy  has  3  degrees  of  freedom.  For  rotational  energy,  the  number  of 
degrees  of  freedom  pT  is  a  function  of  the  number  of  atoms,  na,  given  by, 


na  =  1 
na  —  2 
na  >  3 


The  number  of  effective  vibrational  degrees  of  freedom  can  be  derived  from  the  simple  har¬ 
monic  oscillator  (SHO)  approximation.  Vibrational  degrees  of  freedom  are  therefore  a  function 
of  the  local  temperature,  T,  and  the  species’  characteristic  vibrational  temperature  Qv,  given  by, 

r  2 Qv/T  ua(na  —  1) 

^  ~  e&v/T  _1  2 


7.2  Injection  and  Reflection 

When  a  particle  is  first  injected  into  a  domain  whether  by  initial  conditions,  inflow,  or  surface 
emission,  its  initial  rotational  and  vibrational  energies  must  be  computed.  When  a  particle  hits  an 
accomodating  surface,  its  internal  energies  must  be  recomputed.  These  values  are  sampled  from 
the  equilibrium  distribution  function  for  the  specified  temperature  T.  Internal  energies,  both  ro¬ 
tational,  Er  and  vibrational  Ev ,  as  functions  of  degrees  of  freedom,  pr  and  pv ,  are  computed  as 
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follows.  If  the  number  of  degrees  of  freedom  is  less  than  or  equal  to  two,  the  internal  energy  can 
be  computed  using  a  single  random  number  //, , 


E  =  -  In  (Rl 


ikT 


i<2. 


If  the  number  of  degrees  of  freedom  is  greater  than  two,  internal  energy  is  sampled,  using 
the  acceptance-rejection  method,  from  the  distribution  function, 

1  /  E 


fi.E) 


T(£/2)  Kkl'J 


,-E/kT 


i>2. 


7.3  Energy  Redistribution  During  Collisions 


The  two  mechanisms  for  exchange  of  internal  energy  are  between  translational  and  rotational 
modes  (TR),  and  between  translational,  rotational,  and  vibrational  modes  (TRV).  These  redis¬ 
tributions  take  place  according  to  the  equilibrium  energy  distribution  function  for  a  specified 
number  of  degrees  of  freedom,  given  by, 


f(E)  = 


1 

F^/2) 


\kT  ) 


e~E'kT, 


where  f{E)  is  the  probability  of  the  occurence  of  energy  E,  £  is  the  number  of  degrees  of  free¬ 
dom,  k  is  the  Boltzmann  constant,  and  T  is  the  temperature.  For  an  exchange  between  two  modes 
A  and  B,  with  respective  degrees  of  freedom  (a  and  ,  the  joint  distribution  function  is, 


e-(EA+EB)/kT 

f(EA,EB)  =  r(^/2)  F(^/2)  \kTj 


Ea \uI2~1  (Eb\ 

\kT  ) 


Hb/ 2-1 


If  the  total  energy,  Et  —  EA  +  EB,  is  known, this  can  be  rewritten  as, 

^  ■4’  B>  r((A/2)r  (b  2  y  kT  )  \kTj 

To  sample  from  this  distribution,  the  following  acceptance-rejection  procedure  is  used.  First, 
the  values  xmax  and  fmax  are  computed  as 

1  -  U/ 2 


'Emax  — 


2  -  U! 2  -  tut 2 


f  _/1-r  \W2-lJu/2-l 

J  max  V  -1  ^ max )  ^max 
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A  random  number  R1  is  then  used  to  compute 


/  =  (1  h\)iAl'2  l(//,)Ui/2  ' 


y  =  f/fn 


The  values  /  and  y ,  and  a  random  number  R2  are  recomputed  until  R2  <  y.  Once  this  condition 
has  been  met,  the  energy  Et  is  redistributed  between  the  modes  A  and  B  using. 


Ea  —  R\Et 


Eb  —  Et  —  Ea. 


Once  a  collision  is  selected  to  take  place,  exchanges  may  take  place  between  translational, 
rotataional,  and  vibrational  energy  modes.  A  TR  exchange  takes  place  with  proability  </>#,  and 
a  TRV  exchange  takes  place  with  probability  <f>v.  The  probabilities  (pn  and  <pv  can  be  obtained 
in  two  ways.  They  can  either  be  specified  as  constants  or  computed  as  a  function  of  species 
parameters  and  the  local  translational  temperature. 

In  order  to  compute  the  probabilities  as  functions  of  temperature,  three  additional  parameters 
must  be  stored  for  each  species,  the  limiting  value  of  the  rotational  relaxation  number  ZToo ,  the 
characteristic  temperature  of  the  inter- molecular  potential,  T*,  and  the  effective  excitation  cross 
section  av .  The  TR  exchange  probability  between  species  i  and  j  ,  can  then  be  computed  as 
[Lumpkin91,  Parker59], 


6  \ 
4-2  a) 


'J'-k 

~T 


where  £r  =  +  Ci  '-s  the  number  of  rotational  degrees  of  freedom  in  the  collision,  a  is  the  VHS 

collision  parameter,  and  T  =  T'+TJ  is  the  local  translational  temperature  averaged  between  the 
two  species. 

The  equation  for  the  TRV  exchange  probability,  i f>y ,  was  obtained  from  an  empirical  fit  [Millikan63] 
to  experimental  data  witha  high-temperature  correction  [Park85].  It  can  be  written  as, 


$  = 


e 


tot 


£tot  | 

* ^"re/A /2 


Pc  exp  (AT“1/3  +  B 
kTreJ 


&ref 


(Jv 


2(2  -  a)- 


rnT 


T(2  -  a) 


2(2  -  a) 


r2  kT' 


kTref 


rnT 


kT  \J 7T 


r(2  - 


a 


(2  kT' 

V  mT  J 


1/2  — c 


-1 


mr 
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A  =  FAw.y204v/3 
B  =  FBm3J4Q4J3  -  F0, 

where  Pc  =  101325Pa/Atm  is  a  pres  sure  conversion  factor,  aref  and  Trej  are  the  VHS  reference 
cross  section  and  temperature,  and  mr  is  the  reduced  mass  of  the  collision.  The  constants  FA  — 
2.85  x  lOwK~1kg-1/2,  FB  —  -2.11  x  1015 I\~4Fkg~3/4  and  Fc  —  18.42  are  the  constants 
obtained  by  empirical  fit  [Millikan63] . 


7.3.1  TR  Exchanges 


A  TR  exchange  is  performed  as  follows.  The  total  collision  energy  to  be  redistributed  is  the  sum 
of  translational  and  rotational  energies,  Ec  —  Et  +  Er,  where  Er  =  Ef  +  Ef  is  the  sum  of  the 
rotational  energies  of  the  two  colliding  particles.  The  first  step  in  the  exchange  is  to  distribute  the 
collision  energy  between  translational  and  rotational  energy  modes.  After  checking  the  relative 
velocity  in  a  collision,  the  distribution  function  is  biased,  so  the  number  of  relative  translational 
degrees  of  freedom  must  be  taken  as  it  =  1  —  2a.  The  number  of  rotational  degrees  of  freedom 
is  the  sum  of  rotational  degrees  of  freedom  for  the  two  colliding  particles,  =  if  +  if.  The 
total  collision  energy  is  therefore  first  redistributed  between  £t  and  ir ,  as  described  in  Section 
7.3,  yielding  Et  and  Er ,  respectively. 

The  rotational  energy  is  then  redistributed  between  the  two  particles  with  rotational  degrees 
of  freedom  if  and  if ,  as  described  in  Section  7.3,  yielding  Ef  and  Ef ,  respectively.  The  post¬ 
collision  relative  velocity  g'  is  then  computed  from  the  post-collision  translational  energy  Et 
using  the  reduced  mass  of  the  collision,  rnT, 


/ 

9  = 


7.3.2  TRV  Exchanges 

In  a  TRV  exchange,  the  total  collision  energy  is  Ec  —  Et  +  Ei,  where  Et  is  the  translational 
energy  and  Ei  is  the  internal  energy.  The  internal  energy  is  in  the  sum  of  the  internal  energies 
of  the  two  particles,  Ei  —  Ef  +  Ef .  The  internal  energy  of  a  particle  is  the  sum  of  rotational 
and  vibrational  enegies,  Ef  =  Ef  +  Ef.  Similarly,  the  number  of  degreees  of  freedom  in  the 
various  energy  modes  are,  & ,  £f ,  if ,  if ,  if ,  if ,  and  if . 
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The  total  collision  energy,  Ec,  is  first  distributed  between  translational  and  internal  modes, 
using  £t  and  &,  to  obtain  Et  and  Et.  The  translational  energy  Et  is  used  to  compute  the  post¬ 
collision  relative  velocity,  as  described  above.  The  internal  energy  Et  is  distributed  between 
the  two  particles,  using  degrees  of  freedom  (f  and  (f ,  to  obtain  Ef  and  Ef,  respectively.  The 
internal  energy  of  particle  A,  Ef,  is  then  distributed  between  rotational  and  vibrational  modes, 
with  degrees  of  freedom  (f  and  £f,  to  obtain  rotational  and  vibrational  energies,  Ef  and  Ef, 
respectively.  The  internal  energy  for  particle  B  is  similarly  distributed  between  Ef  and  Ef . 


7.4  Temperature  Calculation 

The  translational  temperature  Tt  for  a  species  is  calcualted  with, 

m(v 2  —  v 2) 

4  =  3 k  ' 

The  rotational  temperature  Tr  is  similarly  computed  from  the  average  rotational  energy  Er,  us¬ 
ing, 

_  ‘2Er 
T  ~  Hr  ' 

The  vibrational  temperature  Tv  is  computed  using  the  characteristic  vibrational  temperature  0 
and  the  average  vibrational  energy  Ev , 

The  total  temperature  T  can  then  be  computed  using, 

rjr  _  tt'Tt  +  irTr  +  £VTV 

6  +  6  +  £v 

where  ^  =  3  is  the  number  of  translational  degrees  of  freedom. 


7.5  Box  Tests 

A  series  of  tests  was  performed  to  validate  this  implementatoin  of  the  Larsen-Borgnakke  en¬ 
ergy  exchange  model.  These  were  designed  to  verify  correct  equilibrium  conditions,  as  well 
as  translational-rotational  (TR)  and  translational-rotational- vibrational  (TRV)  relaxation  rates. 
These  tests  were  performed  for  a  uniform  box  with  specular  walls. 
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Time  (s) 


Figure  15:  TR  Relaxation 


The  first  test  was  for  N2  at  equilibrium  in  a  box,  with  translational,  rotational,  and  vibrational 
temperatures  equal  to  10,000K,  density  1.83e  19  particles  per  cubic  meter,  constant  TR  exchange 
probability  (pu  —  0.2,  and  constant  TRV  exchange  probability  cpv  —  0.02.  The  system  stayed  in 
equilibrium  for  a  =  0,  0.5  and  0.24. 

The  next  case  considered  relaxation  between  translational  and  rotational  modes  for  a  =  0.24. 
Initially,  the  translational  temperature  was  1000K,  the  rotational  temperature  OK,  and  the  vibra¬ 
tional  temperature  OK.  The  TR  exchange  probability  was  constant  at  <f>R  —  0.2,  and  no  TRV 
exchanges  took  place  (<f>v  —  0).  The  results  for  Hawkand  SMILE,  agree  as  shown  in  Figure  15. 

The  next  case  considered  relaxation  between  translational,  rotational,  and  vibrational  modes, 


wit 


l  parameters  as  shown  in  the  table  below. 


Species 

N2 

a 

0.24 

P 

0.735 

<t>R 

0.2 

<t>V 

0.02 

ev 

3390  K 

K 

Trej 

273  K 

K 

U ref 

5.3068  x  10“19 

m2 

Tt 

10000 

K 

Tr 

5000 

K 

Tv 

0 

K 

The  next  case  was  a  mixture  of  Nitrogen  and  Oxygen,  each  with  internal  degrees  of  freedom, 
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Time  (s) 


Figure  16:  TRV  Relaxation 


and  with  TR  and  TRV  exchange  probabilities  as  specified  in  the  table  below. 


Species 

N2-N2 

N2-02 

02-02 

a 

0.24 

0.255 

0.27 

P 

0.735 

0.725 

0.7143 

<f>R 

0.2 

0.2 

0.2 

<f>V 

0.02 

0.02 

0.02 

ev 

3390 

2256 

K 

Trej 

273 

273 

273 

K 

a ref 

5.3068  x  10“19 

5.3069  x  10“19 

5.307  x  10-19 

m2 

Fraction 

.5 

.5 

Tt 

10000 

15000 

K 

Tr 

5000 

7500 

K 

Tv 

0 

0 

K 

As  shown  in  Figure  17,  agreement  was  excellent  between  Hawkand  SMILE. 

The  next  test  considered  variable  exchange  probabilities,  4>r{T)  and  <f>v(T).  As  shown  in 
Figure  18,  the  results  agreed  very  well  with  SMILE. 

7.6  Cylinder  Test 

The  final  Larsen-Borgnakke  test  was  for  supersonic  flow  around  the  cylinder.  Figure  19  shows 
a  comparison  of  temperature  profiles  between  Hawkand  SMILE.  Figures  20  and  21  show  trans¬ 
lational  temperature  flowfields  for  the  two  codes.  Results  agree  to  within  statistical  scatter. 
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Figure  17:  TRV  Mixture 


Figure  18:  Variable  <t>R,<t>v 
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Tvib  (SMILE-OD) 
-•-Trot  (SMILE-OD) 
-©-  Ttrn  (SMILE-OD) 

•  Tvib  (Hawk) 

•  Trot  (Hawk) 

•  Ttrn  (Hawk) _ 


0 . 0 

0.0  0.0  0.0 


Figure  19:  TRV  Cylinder 


r\ 


Figure  20:  Ha wkT ran s  1  at i o n a  1  Temperature 


Figure  21:  SMILE  Translational  Temperature 
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8  General  Collision  Algorithm 


This  section  describes  the  algorithm  for  computing  all  of  the  collisions  in  a  given  cell.  The  vari¬ 
able  recompute  is  used  to  detect  when  the  collision  frequencies  must  be  recomputed.  This  is 
necessary  during  the  first  iteration  of  the  loop,  as  well  as  any  time  a  reaction  takes  place.  The 
total  collision  frequency  v  is  the  sum  of  the  collision  frequencies  between  all  pairs  of  species 
i  and  j . 


1.  recompute=True 

2.  t  =  0 

3.  while  t  <  timestep 

(a)  if  recompute 

i.  let  v  =  J2ij  Vij 

ii.  recompute  =  False 

(b)  T  =  iMM 

(c)  t  =  t  T  r 

(d)  if t  <  timestep 

i.  Select  species  i,  j 

ii.  Select  particles  a,  b 

iii.  if  7;(af)3a6  > 

A.  Perform  Collision 

B.  If  reaction,  recompute  =  True 


The  collision  frequency  between  species  i  and  j  is  computed  using 

—  )  2V  ~  J 

^  _  7 

where  is  the  maximum  value  of  ol3  ( gij )  that  has  been  observed  between  any  pair  of 

colliding  particles  i  and  j,  Jr)\  is  the  ratio  of  real  to  simulated  particles,  and  V  is  the  volume  of 
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the  cell.  Note  that  a  separate  value  of  [<7ijg]max  is  kept  for  each  pair  of  species.  The  species  i 
and  j  are  chosen  with  probability  -A  The  particles  a  and  b  are  selected  randomly  from  the  lists 
for  species  i  and  j . 

Once  two  particles  have  been  selected  for  a  collision,  the  collision  takes  place  according  to 
the  following  algorithm.  A  collision  may  be  reacting  or  non-reacting.  Each  possible  reacting 
collision,  or  reaction,  k,  has  an  associated  cross  section,  <Tk(g),  and  will  occur  with  probability 
.  The  sum  of  cr*  over  all  reactions  is  less  than  or  equal  to  the  total  cross  section  of  the  col¬ 
lision,  at.  Not  all  collisions,  therefore,  are  reacting.  If  a  reaction  takes  place,  internal  energy  is 
exchanged  as  described  in  Section  7.3. 

1.  Select  reaction  k  with  probability 

2.  If  reaction  k  selected 

(a)  Perform  reaction  k 

(b)  It  will  be  necessary  to  recompute  v  and  vij, 

3.  Else  if  no  reaction  selected 

(a)  Exchange  internal  energy  if  necessary. 

(b)  Perform  collision 

4.  Update  particle  velocities  and  species 

If  no  reaction  is  selected,  an  elastic  or  inelastic  collision  takes  place.  With  probability  Ptr,  a 
TR  exchange  will  take  place.  With  probability  Ptrv,  a  TRV  exchange  will  take  place.  If  neither 
occurs,  a  elastic  collision  will  be  performed.  If  a  reaction  is  selected,  the  energy  of  the  reaction 
Ej k  is  used  to  compute  the  post-collision  energy  E'c  from  the  pre-collision  energy  Ec,  E'c  =  Ec  — 
Ek-  The  particles  are  then  moved  to  the  new  species  lists,  and  the  post  collision  velocities  are 
computed. 
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9  Conclusion 


The  H aw  kDSMC  implementation  now  contains  several  collision  models,  Hard  Sphere  (HS), 
Variable  Hard  Sphere  (VHS),  Variable  Soft  Sphere  (VSS),  and  Larsen-Borgnakke.  The  motiva¬ 
tion,  implementation,  and  validation  of  each  of  these  models  has  been  presented.  Results  show 
that  the  implementation  is  correct,  in  comparison  with  another  DSMC  implementation. 

The  tools  described  in  this  paper  are  currently  in  use  by  Intel  and  Tegal  Corporations  for 
design  and  evaluation  of  plasma  reactors,  and  by  the  Institute  for  Defense  Analyses  for  reentry 
calculations.  Future  work  will  include  implementation  and  validation  of  more  sophisticated  re¬ 
actions,  such  as  dissociation  and  recombination.  Reactions  on  surfaces  will  also  be  addressed, 
once  adequate  models  are  available.  Once  these  additional  features  are  in  place,  reliable,  accu¬ 
rate,  and  detailed  reactor  simulations  will  be  possible.  These  in  turn  will  facilitate  the  design  of 
the  next  generation  of  plasma  reactors. 
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